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Abstract. This paper develops a discontinuous Galerkin (DG) finite element 
differential calculus theory for approximating weak derivatives of Sobolev func- 
tions and piecewise Sobolev functions. By introducing numerical one-sided 
derivatives as building blocks, various first and second order numerical op- 
erators such as the gradient, divergence, Hessian, and Laplacian operator are 
defined, and their corresponding calculus rules are established. Among the cal- 
culus rules are product and chain rules, integration by parts formulas and the 
divergence theorem. Approximation properties and the relationship between 
the proposed DG finite element numerical derivatives and some well-known 
finite difference numerical derivative formulas on Cartesian grids are also es- 
tablished. Efficient implementation of the DG finite element numerical differ- 
ential operators is also proposed. Besides independent interest in numerical 
differentiation, the primary motivation and goal of developing the DG finite el- 
ement differential calculus is to solve partial differential equations. It is shown 
that several existing finite element, finite difference and DG methods can be 
rewritten compactly using the proposed DG finite element differential calculus 
framework. Moreover, new DG methods for linear and nonlinear PDEs are 
also obtained from the framework. 



1. Introduction 

Numerical differentiation is an old but basic topic in numerical mathematics. 
Compared to the large amount of literature on numerical integration, numerical 
differentiation is a much less studied topic. Given a differentiable function, the 
available numerical methods for computing its derivatives are indeed very limited. 
There are essentially only two such methods (cf. [Slj)- One method is to approxi- 
mate derivatives by difference quotients. The other is to first approximate the given 
function (or its values at a set of points) by a more simple function (e.g., polyno- 
mial, rational function and piecewise polynomial) and then to use the derivative of 
the approximate function as an approximation to the sought-after derivative. The 
two types of classical methods work well if the given function is sufhciently smooth. 
However, the two classical methods produce large errors or divergent approxima- 
tions if the given function is rough, which is often the case when the function is a 
solution of a linear or nonlinear partial differential equation (PDE). 
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For boundary value and initial-boundary value problems, classical solutions of- 
ten do not exist. Consequently, one has to deal with generalized or weak solutions, 
which are defined using a variational setting for linear and quasilinear PDEs. Al- 
though numerical methods for PDEs implicitly give rise to methods for approxi- 
mating weak derivatives (in fact, combinations of weak derivatives) of the solution 
functions (cf. [31 |6l [HI [32] ) , to the best of our knowledge, there is no systematic 
study and theory in the literature on how to approximate weak derivatives of a 
given (not-so-smooth) function. Moreover, for linear second order PDEs of non- 
divergence form and fully nonlinear PDEs, it is not possible to derive variational 
weak formulations using integration by parts. As a result, weak solution concepts 
for those types of PDEs are different. The best known and most successful one is 
the viscosity solution concept (cf. [131 116] and the references therein). To directly 
approximate viscosity solutions, which in general are only continuous functions, 
one must approximate their derivatives in some appropriately defined sense offline 
(cf. [T71I21]), and then substitute the numerical derivatives for the (formal) deriva- 
tives appearing in the PDEs. Clearly, to make such an intuitive approach work, 
the key is to construct "correct" numerical derivatives and to use them judiciously 
to build numerical schemes. 

This paper addresses the above two fundamental issues. The specific goals of 
this paper are twofold. First, we systematically develop a computational frame- 
work for approximating weak derivatives and a new discontinuous Galerkin (DG) 
finite clement differential calculus theory. Keeping in mind the approximation of 
fully nonlinear PDEs, we introduce locally defined, one-sided numerical derivatives 
for piecewise weakly differentiable functions. Using the newly defined one-sided 
numerical derivatives as building blocks, we then define a host of first and second 
order sided numerical differential operators including the gradient, divergence, curl, 
Hessian, and Laplace operators. To ensure the usefulness and consistency of these 
numerical operators, we establish basic calculus rules for them. Among the rules 
are the product and the chain rule, integration by parts formulas and the diver- 
gence theorem. We establish some approximation properties of the proposed DG 
finite clement numerical derivatives and show that they coincide with well-known 
finite difference derivative formulas on Cartesian grids. Consequently, our DG fi- 
nite element numerical derivatives are natural generalizations of well-known finite 
difference numerical derivatives on general meshes. These results are of indepen- 
dent interest in numerical differentiation. Second, we present some applications of 
the proposed DG finite element differential calculus to build numerical methods for 
linear and nonlinear partial differential equations. This is done based on a very 
simple idea; that is, we replace the (formal) differential operators in the given PDE 
by their corresponding DG finite element numerical operators and project (in the 
sense) the resulting equation onto the DG finite element space Vj^. We show 
that the resulting numerical methods not only recover several existing finite dif- 
ference, finite element and DG methods, but also give rise to some new numerical 
schemes for both linear and nonlinear PDE problems. 

The remainder of this paper is organized as follows. In Section [2] we introduce 
the mesh and space notation used throughout the paper. In Section [3] we give 
the definitions of our DG finite element numerical derivatives and various first and 
second order numerical differential operators. In Section |4| we establish an approx- 
imation property and various calculus rules for the DG finite element numerical 
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derivatives and operators. In Section [5] we discuss the implementation aspects of 
the numerical derivatives and operators. Finally, in Section [6] we present several 
applications of the proposed DG finite element differential calculus to numerical so- 
lutions of prototypical linear and nonlinear PDEs including the Poisson equation, 
the biharmonic equation, the p-Laplace equation, second order linear elliptic PDEs 
in non-divergence form, first order fully nonlinear Hamilton-Jacobi equations, and 
second order fully nonlinear Monge- Ampere equations. 

2. Preliminaries 

Let d be a positive integer, V, C R"* be a bounded open domain, and Th denote 
a locally quasi-uniform and shape-regular partition of r2 . Let denote the set 
of all interior faces/edges oi Th, denote the set of all boundary faces/edges of 
Th, andSh :=£[U£^. 

Let p g [l,oo] and m > be an integer. Define the following piecewise W™'^ 
and piecewise C" spaces with respect to the mesh Th- 

W'^'PiTh) := n ^"'^''iK), C"\Th) J] (^"W- 
Ken Ken 

When p = 2, we set H^^{Th) '■= W"^''^{Th)- We also define the analogous piece- 
wise vector-valued spaces as H'^iT) := [H-^{Th)Y, W'^^PiJ'h) = [W'^''P{Th)Y, 
C^iTh) = [C'^iTh)]'^, and the matrix-valued spaces H'^iTh) := [H''' (Th)]'^'"^ , 
W"''P{Th) [W"''P{Th)]'^'"^, and C"'{Th) = [C"' {Th)]'^'"^ . The piecewise L^-inner 
product over the mesh Th is given by 

{v,w)n ■= ^ vwdx, 
Ken •'^ 

and for a set Sh Eh, the piecewise i^-inncr product over Sh is given by 

eeSn-'" 

Angled brackets without subscripts (•,•) represent the dual pairing between some 
Banach space and its dual. 

For a fixed integer r > 0, we define the standard discontinuous Galerkin (DG) 
finite element space V^^ C IV^'PiTh) C L^ifl) by 

- n ^riK), 

Ken 

where ¥r{K) denotes the set of all polynomials on K with degree not exceeding r. 
The analogous vector-valued and matrix valued DG spaces are given by Vj^ :— [Vj^]''' 
and Vj,^ := [Vj^]'^^''-. In addition, we define 

Vh ■.= w'^\Th)nC'>{Th), 

Vh := [Vh]'^, and Vh := [Vh]'^'"^- We note that Vj^ C Vh- We denote by : 
L^{n) Vj} the projection operator onto , Pj? : [L'^{^)X^ ^ the 
projection operator onto V^'' and : [L-^(ri)] ^ —t^ 1^'' the projection operator 
onto y/'. 
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Let K,K' & Th and e = dK n dK' . Without loss of generality, we assume that 
the global labeling number of K is smaller than that of K' . We then introduce the 
following standard jump and average notations across the face/edge e: 

[v] := v\k — v\k' on e e f/J, [t;] := u on e € f^, 



{'"}'■= ^{""Ik + v\k') oneGSl, {v} := v on e G £, 



for V e Vh- We also define rig := nK\e = —nK'\e as the unit normal on e. 

3. Definitions of discrete differential operators 

Let V e Vh- For e e f^, that is, e = dK n dK' e £^ for some K, K' e %, we 
write rie = [ni^\ nf \ . . . ,ni'*^)* to be the unit normal of e. We then define the 
following three trace operators on e in the direction Xi: 

(3.1) Qi{v){x):- 



(3.2) 
(3.3) 



Qtiv){x) 



lim 


^^(y) 


if 




<o, 


lim 


v{y) 


if 




>o, 


lim 

y€K> 

y— 


v{y) 


if 




<o, 


lim 

yeif 


v{y) 


if 




>o, 



-^[Q-;{v){x) + Qt{v){x) 



for any .x G e and z = l,2,...,(i. Wc note that and can be regarded 
respectively as the "left" and "right" limit of v at a; G e in the direction of Xi. If 
e e f/f , we simply let 

(3.4) Qi{v){x) = Qt{v){x) = Qi{v){x) := lim v{y) Va; e e. 



yen 

y— >x 



RemEirk 3.1. On an interior edge e G E^, we may alternatively write 



(3.5) Qf{v) = {v}±-sgn{n'i^)[v], where sgn(nW) 



ifn^e ^ > 0, 
ifn^e^ < 0. 



With the help of the trace operators , and Qi we are ready to introduce 
our discrete partial derivative operators d'^^,,d^^,, dh,xi ■ V/i Vj^. 



Definition 3.1. For any v G Vh, we define the discrete partial derivatives 
K-,^' ^h,x,v e V,^ by 

(3.6) {d^^^^v,^h)r, ■■= {Qfivh^'\ - («,5,.^/.)r. 



h,Xi 



(3.7) 



for i = 1, 2, . . . , d. Here, dx^ denotes the usual (weak) partial derivative operator in 
the direction Xi, n^*-* is the piecewise constant function satisfying and 
7~ and 7^ are piecewise constants with respect to the set of interior edges. 
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In addition, we define the discrete partial derivatives when boundary data is 
provided. 

Definition 3.2. Let g G L^{drt) be given. Then for any v G Vh, we define the 
discrete partial derivatives 9^ ^ d^'f^v, 9^ ^. f G by 

(3.8) {d^'lv, <^,)^^ := {dl^v, ip„)r, + ((.g - v)n(^\^,,)^^ e K^ 

Remarks 3.1. 

(a) Since every function v ^ Vh has a well-defined trace in L^{dK) and every 
function ip^ S ^'^^ well-defined trace in L°°{dK) for all K G Th, the 
last term on the right-hand side of (3.6) is well defined. 

(b) Since Vj} is a totally discontinuous piecewise polynomial space, the discrete 
derivatives d^^ v can also be written in their equivalent local versions: 



K 



eCdK\dn 

for i = 1,2, ... ,d and K eTh- Here, ^f^^ = lf\e- 

(c) The discrete derivatives dj^ ^ ''^'^h '^"'^ (^h,xiV can be regarded, respec- 
tively, as "left", "right", and "central" discrete partial derivatives of v with 
respect to Xi. The definitions are analogous to the weak derivative defini- 
tion. 

(d) We note that the discrete one-sided partial derivatives are defined for func- 
tions in the space Vh, in particular, for functions in the DG finite element 
space T/'' C Vh- 

(e) By the identity (3.5), we have 

{dl^v,Vh)^,^ = {{vW\ [vh])^ + Ulf ± \\n^'^)[v], [Vh])^, - {v,d,^Vh)^^. 

Therefore, d^^v — df^,^,v provided — = — for all e G 
Definition 3.3. We define the following first order discrete operators: 
(3-lla) ^hV-^{dt,v,d^^^^v,--- 

(3.11b) V^z;:=i(v^t; + V+i>), 

(3-llc) div±t; d^^^^v, + d^^^^v^ + ■■■+ d^^v^, 

(3. lid) div/jt) := ^(div^v + div^i;) 

for any v G V/i and v = (wi,W2, • • • ,i'd)* G V/j. . In addition, we define the discrete 
curl operators for v G Vh o,nd v G V/,.; 

(3.11e) curl^t; := (^^,^^3 - dl^,V2,dj^^^,^vi - d^^^v^,d^,^^V2 - ^L.^i^ 
when d ~ i, and 

(3.111) curl^t; :=9±^^V2-5±^^vi, cuy\^v:= (d^^^^v,-d^^^^v^ 
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when d — 2. We also set 

(S.llg) curl/i?; - ^curl^t; + curl^t;^ , curl hV = - (curl ^t; + curl ^v) . 

When boundary data g G L^{dn) is provided, the analogous operators are given by 
(3.12a) V±^«:= (C'"' ' ' ' ' 

(3.12b) V;,,,«:=^(v^'«,; + V^^i; 

In addition, for given g — {gi,g2, ■ ■ ■ ,gdY G L^{dVt), we set 
(3.12c) div±^^ ^^:^lv^ + d^^v^ + • • • + 

(3.12d) ^\-^h,gV := ^(div^gt; + div^gv). 

Definition 3.4. Denote by and Df^ the transposes of the operators and 
V/i, respectively; that is, D^v — (Vjw)* and D^v = (V/iu)*. We then define the 
following second order discrete operators: 

(3.13a) D-^v := D-Vf^v, D+^v D+V^v, 

(3.13b) A,;±z;:=tr(i?^±«), A+±« tr(i^+±t;), 

(3.13c) Dlv.^DhWhV, Ahv:= ti-{Dlv) 

for any w G V^. Here, tr(-) denotes the matrix trace operator. When boundary data 
g G L^{dfl) is given, we define 

(3.14a) ^,;Vt^, D+^^ := D+^l^v, 

(3-14b) tr(l?-», A+±^; tr(i?+>), 

(3.14c) i^^^gW DhVh,gV, Ah^gV tr(i?,2,_gw). 

Remarks 3.2. 

(a) The matrix-valued functions Dj^~v,D'^^v,Dj^~^v and D'^~v define four 
copies of discrete d X d Hessian matrices of the function v. Similarly, 
A^~w, A^^u, A^"*"?; and A~^~v are four copies of discrete Laplacians of v. 

(b) Since V^w, V^w G Vj^ C V/i, the above second order discrete operators 
are well-defined. 

(c) It is simp le to see that the definitions of the discrete Laplacians given in 



Definition 3.4 are equivalent to Aj^^v = div^ V^w, A^^w = div^Vjw and 



AhV = divhVhV. 



4. Properties of discrete differential operators 

In this section we shall establish a number of properties for the DG finite element 
(FE) discrete derivatives defined in Definitions 3.1 and 3.3 We start with some 
characterization results for the DG FE discrete derivatives. We then establish an 
approximation property followed by several basic calculus rules such as product 
and chain rules, the integration by parts formula, the divergence theorem and 
the relationship of the DG FE discrete derivatives with standard finite difference 
derivative formulas. 
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4.1. Characterization of DG finite element derivatives. For any v e H^iVl), 
we have Qf{v)\e = v\e and [w] = on e e Hence, 



h 



Integration by parts element-wise immediately yields 

(4.1) V,^)^^ = (^;,)^^ V(^;, G K^ 

and therefore 

(4.2) fe.^,</'/.)K = y'Phe¥r{K),yKeTh. 

Hence, we have the following proposition. 

Proposition 4.1. For any v e V;i n H^{il,), d^^ v and dh.xi^ coincide with the 

-projection of dx^v onto Vj.^ . We write d^^,v — dh,xiV = Vl^dxiV, where V!^ 
denotes the projection onto Vj} . 

From the above proposition, we easily derive the following corollary. 

Corollary 4.1. Let v,, G V^''nH^{n) with 0<£<r + l. Then dJ^.^^Vh = dj^^^^Vh = 
dh,XiVh = dx^Vh € Ve-i' where V!^i :— {Q}, the set with only the zero function. 

For an arbitrary piecewise function v e Vh, the above characterization results 
are not expected to hold in general because dx^v may not exist. Below we shall 
derive a similar characterization to (4.1 ) for d^^ v and dh.xiV when v is an arbitrary 
function in Vh- 

For any v &Vh, let T^x^v denote the distributional derivative of v with respect to 
Xi (which does exist). Let 5 C SI be a (d— l)-dimensional continuous and bounded 
surface. We define the delta function (5(E!, g, x) of variable strength supported on S 
by (cf. M) 

(4.3) {6{E,g,-),^):^ [ g{sMx{s))ds V(p G C°(r!), 



where x{s) G ^. We also extend the above definition to test fvmctions from Vj} as 
follows 

(4.4) {S{E,g,-),^h):= l^g{s){M<s))}ds V^, G V;^ 

Using 6{'E.,g,x) we give the following characterization for Vx^v. 
Lemma 4.1. For every u G V/i there holds the following representation formula: 

(4.5) Vx,v = ^ dx^vxK - ^ n'^e^6{e,[v],x) fora.e.xefl, 

Ken ee£^ 

where xk denotes the characteristic function supported on K . 
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Proof. By the definition of Vx^v and integration by parts on each if e fi, we have 
for any ip e (n) n (fl) 

K^Th 

^ n«(S(e,H,-),<^) + ( ^ d.j:,vxK.V 
eeei Ken 

Here we have used the fact that v E W^'^{K) for every K e Th- Clearly the above 
identity infers (4.5). The proof is complete. □ 



Proposition 4.2. Set 7 — m (3.6). For any v € Vh, dh^xi^ coincides with the 
-projection ofUx^v onto Vj} in the sense that 

(4.6) {dh,x,v,iph)j-^^{'Dx,v,iph) yifih^V, 



h 

r 5 



where the right hand-side is understood according to (4.4). We write dh,xiV 



V^Vx^v. 



Proof. By the definition of dh^xi^ and using the integration by parts formula for 
piecewise functions, we get, for i = 1, 2, . . . , d. 



e(^£{ Ken 



Here we have used (4.4) and (4.5) to get the final two equalities. The proof is 
complete. □ 

Remark 4.1. Define the "lifting operator" Ch/i ■ L^{£h) -> Vj^ by 
(4.7) (/:M«,^/Or. (E "^'^-^(^'Hr),^/.) = 



eeei 



Ken 



Then we have 
(4.8) 



Finally, combining Propositions |4.1| and |4.2| and the well-known limiting charac- 
terization theorem of distributional derivatives (cf. [251 Theorem 6.32]), we obtain 
another characterization for our DG FE derivatives. 

Proposition 4.3. Set 7± = m (|3^. Then for any v Cz Vh, there exists a 
sequence of functions {vj}j>i C C^{fl) such that, for i = 1, 2, • • • , d, 

(i) Vj V as j ^ CO in L^(fl). 
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(ii) dh^xiVj — > V^'Dx^v as j —> oo weakly in in the sense that 
(4.9) \iin{dh,x^vj,^h)r ^{K"Dx^v,^h)r "i^h^V^. 

j—^OQ ' I h ^ ' I h 



Proof. Let {pj}j>i denote a symmetric moUifier (or approximate identity) with 
compact support. For any d € Vh C L^{^), set Vj := v * pj € C^{Q), the 
convolution of v and pj. Then it is easy to check that the sequence {vj} fulfills the 
properties (i) and (ii). The proof is complete. □ 



4.2. Approximation properties. By Proposition |4. 1| we readily have the follow- 
ing approximation properties for the DG FE derivatives. 

Theorem 4.1. For any w G V/i n H^{fl) there holds the following inequalities: 

(4.10) 115,,^; - OIxAlHU) < Wd.^v - Ox^Ml^u) ^i'h e K+i, 

(4.11) \\dxiV - dl xM\L'^iK) < \\dx,V - dx,^h\\L^{K) ^Iph e Pr+1, 

where * can take +, — and empty value. 
Proof. Notice that ( |4.1[ ) can be rewritten as 

{dx^v ~ dlx.v, ^h)r,^ =0 y^he v^. 

Hence, for any iph € Vj^+i, there holds 

\\dx,v - dl^xM\h{U) = (^^.^ ~ ^h.x,^, dx,v - d'lx^v)^^ 
= {dx^v - dl x,v,dx,v).^^^ 
= {dx,v - dl x,'^^ dx,v - dxii>h)j-^ 
< \\dx,v ~ dx^v\\L2(r^)\\dx,v - dx^iphWL^u)- 



The identity (4.10) follows from the above inequality because dx^Vl^i C . The 
second inequality (4.11) is obtained by similar arguments. □ 

It is also possible to derive some error estimates for Vx-V — 9,* ^.u with general 
V £ Vh in some weaker (than L^) norm. However, since the derivation is lengthy 
and technical, we leave it to the interested reader to exploit. 



4.3. Product rule and chain rule. The product rule and the chain rule are 
two very basic properties of classical derivatives and weak derivatives (cf. \2.1\ 
Chapter 7]). The goal of this subsection is to establish both of the rules for our DG 
FE derivatives. As expected, these discrete rules take different forms from their 
continuous counterparts. 

First, we consider the case when functions are from the space Vh H C''(fi). In 
this case, we have the following product rule and chain rule. 

Theorem 4.2. Let F € C^(R),F' e L°^(R). For u,v eVkH C°{n), there holds, 
for i ^ 1,2, - ■ ■ ,d, 

(4.12) dh,xAuv)= {udx, V + vdx, u) , 

(4.13) dh,xj{u) = V'^{F'{u)dx^u). 
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Proof. Notice that uv eVhH C°(fi) for any u,v e Vh n C^{Q). By Proposition 
and the product rule for weak derivatives (cf. [21, Chapter 7]) we get 



4.1 



Hence, ( |4.12[ ) holds by the definition of "P^. 

Similarly, since F{u) G V,, n C"(f7) for u e V/^ n C"{n) we have 

which then infers ( 4.12^ . The proof is complete. 



□ 



An immediate consequence of Theorem 4.2 and Corollary 4.1 is the following 
corollary. 

Corollary 4.2. Let F e C^{Tl),F' e L°°{R). For u,v € n C"(f7) with < 
i — 1 < r, there holds, for i ~ 1,2, ■ ■ ■ ,d, 



(4.14) 
(4.15) 



dh^x^uv) = 'P'^{udh,x,v + vdh,xiu), 
dh,x,F{u)^V'^{F'{u)dn,x^u). 



Next, we consider the general case when functions are from the space Vh- As 
expected, the discrete product and chain rules appear in more complicated forms. 

Theorem 4.3. Let F e C^R) and F' e i°°(R). For u,v E Vh, there holds for 
i = 1, 2, • • • , d, 



(4.16) 
(4.17) 



dh^x,F{u)=V^{F^V^^u), 



where T^x^v is given by (4.5) and v represents a modification of v which assumes 
the same value in each element K G 77i and takes the average value of v on each 
interior edge e S Thus, 

(4.18) uVx^v — udxiVXK — n^f^ 5{e,{u}[v],x) for a.e. a: e fi. 

Proof. Without loss of the generality, we assume F g C°°(R). For u,v e Vh, let 
{vj}, {uj} C C§°{n,) be defined in Proposition 4.3 Then, by Theorem 4.2 we have 



(4.19) 
(4.20) 



dh.xiiujVj) = VriujdxiVj + VjdxiUj), 
dh,x,F{uj)^V':{F'{uj)dx,Uj). 



It follows from ( [4.9^ and (4.6) that for any iph G , 

(4.21) \\Ta{dh,xSu,Vj),iph)r = {V^Vx^{uv),^h)r 

= {Vxi{uv),iph) = {dh,xiiuv),(ph)^^, 

(4.22) lim {dh.x,F{u^),iph)^^ - {V'^V^^F{u),iph) ^ 

= {'DxiF{u),iph) = {dh,xiF{u),(ph)j-^ 
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On the other hand, we have 



(4.23) hm (V'rUjdx^Vj + Vjdx^Uj) , iph) ^ 

j—^OO ' ' ' h 



(4.24) 



hm (ujdx^Vj + Vjdx^Uj,iph)^ 
\im ('Pr(F'{uj)dx,Uj),iph)^ = hm (F'{uj)dxiUj,iph)^ 



= {F'{u) V,M,^h)- 

Then, ( |4.16p follows f rom c ombi ning j4.19| ), ( |4.21| and ( |4.23[ ), and ( |4.17| ) follows 
from combining (4.20), (4.22| and (4.24). The proof is complete. □ 



Remark 4.2. T/ie smoothness assumption on F in Theorem 4.3 may be weak- 
ened to F £ C'^'^(R) by following the techniques used in pp. We leave the general- 
ization to the interested reader to exploit. 

4.4. Integration by parts formula and discrete divergence theorems. In 

this subsection, we derive integration by parts formulas for the discrete differential 
operators that resemble the standard integration by parts formula in the continuous 
setting. These results will play a crucial role in developing DG methods for a variety 
of PDE problems as well as in relating these methods to other existing DG methods 
in the literature. 

Theorem 4.4. Suppose that = that is, ~ ^l7e ^ ^ ^h- Then 

the integration by parts formulas 



(4.25) 
and 

(4.26) {dh,xiVh,(Ph)j-^ = -{vh,dh,xi(Ph)j-^ + {vh,fh'n'^''^)£B 

hold for all Vh, fh G Vj^- 



Proof. It sufhces to show (4.25), since (4.26) trivially follows from (4.251. 



By ( |3.6| and integration by parts, we obtain 

h 

Using the identity ( |3.5[ ), we have 

(4.27) {dl^^v.Vh)r, = {d..v,Vh)r, ± ^(H, [v'Jsgn(nW)n«)^, 
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Now let V = Vfi ^ Vj\ Then by (3.6), (3.5) and rearranging terms, we have 
(4.28) - ^{dl^^^^,Vh)^^^ + {{^^}, [vhW'^)^, 

Using the identity ( 4.28[ ) in equation ( 4.27| , we obtain 

(4-29) {dj^^^^Vh,Vh)^^ = ~{vh,dl^^^h)r, + {vh.'Phn'^'^)^^ 



The identity (4.25) then easily follows. 



□ 



~li for all i = \,2, . . . d. Then the integration 



Theorem 4.5. Suppose that 7^ 
hy parts formulas 

(4.30) {dbf^Vh^'fh)-^^ = -{vh,V'^iph)^^+ {vh-n,iph)^E 
and 

(4.31) (divhVh,iph)j-^ = -{vh,'^hiPh)j-^ + {vh ■ n,iph)gi 
hold for all Vh G and iph e Vj^ . 



Proof. If 7+ = --f-, then by (3.11c), (4.25) and (3.11a), we have 



(divjvft, iph)j-^ = {^h,x.'"h,i, fh) 



n 



1=1 

d 



i=l 



Formula (4.31) is obtained similarly. 



□ 



Theorem 4.6. The formal adjoint of the operator div^ (resp., div^jj is — V^q 

(resp., — V/i.oj with respect to the inner product (•, provided = —J^ for all 
i ^ 1,2, ... ,d; that is, 

(4-32) {diY'^Vh,iph)r^ = -(t;,„Vjo^ft)r^, 

(4.33) {div hVh,'Ph)n = -{""h,^ hfiVhWh 

for all Vfi £ Vj'', fh G ■ ^'^ addition, if 'y^ — ~7.j then the formal adjoint of 
the operator div Q (reps., div/,,_oJ is — (resp., — V^j. 

Proof. This result immediately follows from Theorem |4 . 5 1 and the identities 

{v, yto^h)r^ = {v, Vh^h)^^ -{v-n, iph)^s, 
(div^Qi;, iph)j-^ = (div^f , iph)^^ - {v ■ n, iph) 
for all V £Vh and iph £ Vj\ □ 
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Theorem 4.7. Suppose that 7+ 

(4.34a) -(A±+t;,^,Or. 
(4.34b) -{^f-v,^H 



Th 



-7 . We then have 



and 
(4.34c) 



for all V £ Vh and iph G Vj}. In addition, under the same hypotheses on jf, we 
have 



(4.35a) 
(4.35b) 

and 

(4.35c) 



n 



{^h,gV,'Vhfi(ph, 



n 



Proof. These formulas foUow from Theorems 

.:9 



Vh = V^gV e Vj' (cf. Remark |3^). 



4.5 



4.6 



with Vh = V^w e Vj^ and 

□ 



4.5. Relationships with finite difference operators on Cartesian grids. We 

now show that when 7^ ~ the operators 9^^. and dh.xi are natural extensions 
(on general grids) of the backward/forward and central difference operators defined 
on Cartesian grids. 

Suppose Th is a, rectangular mesh over Q which is aligned with the underlying 
Cartesian coordinate system. Let hi denote the mesh size of Th in the direction of 
Xi. Notice that when r ~ 0, the function d^^.v is a piecewise constant function 
over the mesh Th- Setting f — 1 in (3.10) we get 



(4.36) 



We only consider the two dimensional case. Then K has four edges with = 
(—1,0)*, (1,0)*, (0,1)*, (0,-1)*. Let be a grid function over Th and Vij denote 
the value of v on the {i,j) cell/element. By (4.36) we obtain 

W, + lj 



(4.37) 
(4.38) 



hi 



h,X2 "^3 



i-lj 



hi 



Consequently, we have 
(4.39) dh^xiVij = 



2hi 



h,X2 



dh'x.Vi 



Vij + 1 — Vij^l 



2/i, 



Hence, d^^, and d^^, coincide respectively with backward and forward difference 
operators in the direction Xi, while dh.xi results in the central difference operator 
in the direction x^. 

Using the operators 9^^. as building blocks, we can also build DC FE operators 
that are natural extensions of the standard finite difference approximations for 
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second order partial derivatives. Frora (4.37) and (4.38), we have 



(4.40) Ci^^",.."»J=^/:,.id^^. 



(4.41) dh^xidh,xiVij = 



Vi~2j - "^Vjj + Vi+2j 



(4.42) 



dtx^dh,x2^i3 + \x^9lx2^^j 



Vi^ij + Vij^i ~ Vi^ij+i ~ 2vij - Ui+ij-i + Uij+i + Ui+ 



(4.43) dh^x,dh,x 



2hih2 



Vi+lj + l — Wj+lj-l — Wi-lj + 1 + Vi-lj-1 
4/li/l2 



Thus, for k ^ £, the discrete differential operator (9^^ '^hx. ~^ x ^hx )/2 coin- 
cides with the standard second order 3-point central difference operator for non- 
mixed second order derivatives, dh,Xk(^h,Xk coincides with the standard second order 
3-point central difference operator with mesh size 2h, {d'^x^^h xe + Xk^h xj/'^ 
coincides with the second order 7-point central difference operator for mixed second 
order derivatives, and dh.x^dh.xi coincides with the standard second order central 
difference operator for mixed second order derivatives. 



5. Implementation Aspects 

In this section we explain how the discrete partial derivatives are computed. 
Denote by K :— {^x £ R'* : > 0, X^iLi -'^l reference simplex, and let 
{'t'h^}%i denote a basis of Vr{k). Here = dimP^ = C'^'^)- For example, we 
could take the basis to be the space of monomials of degree less than or equal to 
r: {x°' : \a\ < r}. For K S Th, let Fk : K ^ K he the affine mapping from 
K onto and and let ^'f^'^^ : A' — > R be defined by ip^l'^\x) = iffp{x), where 
X — Fk{x) £ K. It is then easy to see that is a basis of Vr{K). We 

then define the mass matrix Mk G R"'-><"'- associated with K as 

By a change of variables we easily find Mk — \<lei{DFK)\M — d\\K\M , where 
M is the mass matrix associated with K, DFk is the Jacobian of the mapping Fk 
and \K\ is the d-dimensional volume of the simplex K . 

Next, given v S Vu, write d^,^^v\^ = E"=i "j'^^^Vh^'^^ e Vr{K) with af'^'^ = 
af^K ^ ^ (* — j = l,2,...,nr). We then define the vector — 

bf,,{v) e R"'- by 

+ E 7t(H'bL-"'"])e' J = l,2,...,n,. 

e<zaK\dQ. 
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Then by (3.101, the coefficients {af^Yj^i are uniquely determined by the hnear 

>,± - h± ,„I,„.o ^± - r^±(l) ^±(2) ±("-)^t 



equation AfA-a- ~ b- , where a- = {a^ ^ , . . . . Equivalently, we 

have = gji^M^^b,*, where M'^ 
matrix which can be computed offline. 

Now a V — v/i G Vj}, then we mf 
constants P'-^^ = f3''^-^'> G R. We then define the matrix Qf = Qf^ e R"'->^"'- by 



have at = -^jjjj^M^^bf , where is the inverse matrix of the reference mass 

;h can 

Now if V — Vfi € V^, then we may write Vh\j^ = X]j=i Z^'"* for some 



eCdK\dn 



for TO = 1, 2, . . . , rir- Again, writing d,^ ^.v\^ = J2j=i '^i ' ^^^'^ ^hat 



the coefficients satisfy = ^^M-^Q^f3, where ^ = (/3(i), /3(2)^ . . . ^ /^("r))*. 



6. Applications 

In this section, we apply our DG FE differential calculus framework to construct- 
ing numerical methods for several different types of PDEs. Essentially, we replace 
the (continuous) differential operator by its discrete counterpart. In addition, we 
add a stability term which ensures that the discrete problem is well-posed and also 
enforces the given boundary conditions weakly within the variational formulation. 
Throughout the section, we assume the penalty parameters 7^ (which appear in the 
definition of the discrete differential operators) are zero in the discussion below. 



6.1. Second order linear elliptic PDEs. 



6.1.1. The Poisson equation. As our first example, we consider the simplest 
linear second order PDE, the Poisson equation with Dirichlet boundary conditions: 



(6.1a) -Am = / in f7, 

(6.1b) u = g on Sfi, 

where / € Li^lQ) and g € LP'{d^) are two given functions. We then consider the 
following discrete version of (|6.1[): Find Uh € such that 



(6.2) -^KgUh+jh,g{'^h)^V'^f. 

Here, jh.g{-) ■ '^h ^ Vj^ the unique operator satisfying 

and rji is a penalty parameter that is piecewise constant with respect to the set of 
edges. 

Problem (6.2) has several interpretations. On the one hand, by the definition of 
the discrete Laplace operator, the problem is the twofold saddle point problem 

fh = DhQh, Qh = '^h.gUh, 
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with Qh e V^^ and fh G . Namely, problem (6.2) is equivalent to finding Uh G V/i, 
Qh e Vf^ and fh e such that 

(6.4a) (q,i, t/Jt^^ = ({"/i}, [t/i] ■ - (w^, divT/j)7; + {g, ■ n)^^ Vr,j e F^'', 

(6.4b) {fh,P'h)^^ = {{qh}, [P'h])i;^ - (QhA^'^P'hh^ ^P'h e VJI, 

(6.4c) - (tr(f,,), (^h.)^^ + (j,,,g(u,,.), (p,,)^^ = (/, v)r^ yiph e v;''. 



On the other hand, by Theorem |4. 5 [ we can write problem (6.2) in its primal form: 
find Uh e Vj} satisfying 

(6.5) {^h.^gUh,^h'Ph)TH - {^h,gUh ' n, 'Ph) + {jh,g(uh), 'Ph)j-^^ (/, 'Ph)n 

for all ip £ . Alternatively, by Theorem 4.6 we may write 

{"^ h,gUh,'^ hfi'-PhWh + {jh,g{uh), fh)j-^^ = if, fh)n ^fh S V^^ . 

Finally, in the case 5 = 0, problem (6.5) can be viewed as finding a minimizer of 
the functional 



(6.6) 



^ / I'^hfiVhl'^ dx + y2 [ yyilKll^c^s- [ fvhdx 



over all Vh € V^^ . 

We now discuss the well-posedness of problem (6.2) as well as relate the dis- 
cretization to other DG schemes. Let qh = ^h,gUh, that is, qh € Vj^ satisfies 
(6.4a). Then by (6.5), we have 

(6.7) {qh,^h(ph)n - {qh ■ n,iph)^B + {jh,g{uh),vh)j-^^ = {f,fh)n 

for all iph € Vj^. By the definition of the discrete gradient and integration by parts, 
we have 

{qh,^hfh)TK = {[qh] ■ n, Wh})^^ - (divq,,, Lph, 
= {qh,"^Vh)Th -{{qh}-n, VPh])^,^- 
Using this identity in ( |6.7[ ), we find 

(6.8) {qh,"^^h)n - {{qh} ■ n, [y^h])^^ + {jKgi^h), Vh)^^^ = (/, ^h)Tu Vcp G . 



In summary, problem \Q.2\ is equivalent to the mixed formulation (6.4a), (6.8). 
This formulation is nothing more than the local discontinuous Galerkin (LDG) 
method [TOlEj. We then have (see, e.g., [8]) 

Theorem 6.1. Let r > 1, rji > Q and 7i = (z = 1, 2, . . . , d). Then there exists a 
unique Uh G satisfying (6.2). Moreover, if rji = 0{h^^) and if u G H'^~^'^{n), 
there holds 

II" - UhWL^n) + h\\Wu - Wh,gUh\\L^n) < Ch'^+^\\u\\Hr-+2(^Qy 

Remark 6.1. A similar methodology can be used to construct DG schemes for the 
Neumann problem 

Ou 

(6.9) — Au = / in ft, —— ^ q on dft. 

on 

In this case, the DG method is to find Uh G satisfying 

-diVh.q'^hUh + jh{uh) = V^f, 
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where jh{uh) ■ is the operator satisfying [jh{v),iph) = ('7iH; b/iDfi 

for all (pfi E Vj}, and q — qn. We again recover the LDG method for the Neumann 
problem ([6^ (cf. 

6.1.2. The DWDG method for the Poisson problem. In this subsection we 
formulate a new DG method for the Poisson problem (6.1) that inherits better 
stability than the LDG method described above. The new scheme, called the sym- 
metric dual-wind discontinuous Galerkin (DWDG) method, is simply given by 



.10) 



+ Jh,giuh)^V';f, 



where jh.g{uh) is defined by (6.3). By Theorem 4.7 problem (6.10) is equivalent to 
finding Uh £ Vj} such that 



.11) 



\{i'^h.g'^h,'^tfiVh)n + i'^h,g'^h,'^h,oVh)n) +jh,g{^h) = {f,Vh)n 



for all Vfi G Vj} . Equivalently, in the case g 
unique minimizer of the functional 



0, problem (6.10) asks to find the 



1 



dx 



^ 2 



»7l \Vh\ 



' ds 



fvhdx 



over all Vh £ (compare to (6.6)). A complete convergence analysis of the sym- 
metric DWDG method for the Poisson problem is presented in [55]. Here, we 
summarize the main results, namely, well-posedness and optimal rates of conver- 
gence. 

Theorem 6.2 ([H]). Set 7niin := miUegf^ h~^r]i{e). Suppose that there exists at 
least one simplex K £ Th with exactly one boundary edge/face. Then there exists 
a unique solution to (6.10) provided ^irnn ^ 0. Furthermore, if the triangulation is 
quasi-uniform, and if each simplex in the triangulation has at most one boundary 
face/ edge, then there exists a constant C* > independent of h and rji such that 
problem (6.10) has a unique solution provided jniin > — C*. 



Remark 6.2. We emphasize that problem (6.10) is well-posed without added penalty 
terms. As far as we are aware, this is the first symmetric DG method that has this 
property in any dimension (cf. |27l 1251 15] ). 



Theorem 6.3 ([26 ). Let Uh be the solution to (6.10), u £ H^^^{^) be the solution 
to (6.1) and 7max — maxeg^^ h~^T]i{e). Then Uh satisfies the following estimate 
provided ^i^iin > 0.' 



(6.12) 



/ 1 \ 2 

^ Y 7miTi ^ 



and if the triangulation is quasi-uniform and 7,„i 

(6.13) \\u~Uh\\L-(n) < C/i^+'(Vl7min| + 



> —Cs,, then there holds 

1 n2 



VC* + 7n 



l"lif= + i(0), 



where C denotes a generic positive constant independent ofh, and C* is the positive 
constant from Theorem |6.2[ 
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Remark 6.3. In light of (4.40) and (4.411, we can see that when approximat- 
ing with piecewise constant basis functions on Cartesian grids, the DWDG method 
coincides with the standard finite difference method for Poisson's equation while 
the LDG method coincides with a staggered finite difference method for Poisson's 
equation that uses coarser second derivative approximations. 

6.2. Fourth order linear PDEs. In this subsection we show how to use the dis- 
crete differential calculus to develop DG methods for fourth order linear PDEs. For 
simplicity, we focus our derivation to the model biharmonic problem with clamped 
boundary conditions: 

(6.14a) A^u = f inQ, 

(6.14b) u^g ondn, 

du 

(6.14c) 1^ = 9 on on, 

on 

where we assume that / € L^(il) and g,q £ L^{dft). The biharmonic problem 
with other boundary conditions (e.g., simply supported or Cahn-Hilliard type) are 
briefly discussed at the end of the section. 

To develop DG methods using the discrete differential machinery, we first need to 
define some additional discrete operators corresponding to the biharmonic operator. 
Given the two functions g and q defined on the boundary, we set 

(6.15) Ah,g,q ■ = div,i_qV/i^g with q = qn, 
and 

(6.16) • A/i(A,i,g^q) div/iV/idiv,i,qVh,g. 

In addition, we define the operator r,,,,(-) : W^^^{Th) n C^iTh) V/' by 

{rh,qiv),iph)j-^ = {r]2[dv/dn],[diph/dn])^j^ + {ri2{dv / dn - q),diphl dn) 

which we will use to enforce the Neumann boundary condition weakly in the DG 
formulation. 

The DG method for (6.14) is then defined as seeking Uh e such that 



•17) K,g,q^h + ]h{uh) + rh{uh) = V;f^ 



Similar to the discussion in Section |6.1[ we may write the DG method in various 
mixed forms (with up to six unknown variables). Instead, wc focus mainly on the 
primal formulation. 

By ( |6.15[ ) and Theorem |4.6[ we have 

{AlgqUh,(ph)j-^ = -{^ hAh,g^qUh,V hfliph)^^ = {A.h,g,qUh, dlYhfl'^ h,OVh) j-^ 

Thus, we may write ( |6.17[ ) in its primal formulation as follows: Find Uh & Vj} such 
that 

(6.18) {Ah^g,qUh, Ah,0,0'Ph)j-^ + {jh,g{Uh), fh)-j-^ + {rh,q{Uh),'Ph)j-^ 
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Remark 6.4. The DG method (6.17) closely resembles the local continuous discon- 
tinuous Galerkin (LCDG) method proposed in 22J (also see [331 J. Here, the authors 
consider a mixed formulation of the biharmonic problem with the Hessian of u as 
an additional unknown. The derivation of the LGDG method closely resembles the 
derivation of the LDG method for the Poisson problem; the main difference is that, 
as the name suggests, the LGDG method uses continuous finite element spaces. 

Theorem 6.4. Suppose that 771 > and 772 is non-negative. Then there exists a 
unique G V^'' satisfying (6.171. 

Proof. Since the problem is finite dimensional and linear, it suffices to show that if 
/ = 0, (7 = and q = Q, then the solution is identically zero. 

To this end, we set = ^hfiUh, Vh = divh^qh, and Zh = VhVh so that d\VhZh + 
jh,o{uh) + 'rh,oiuh) = 0. To ease notation, we define the bilinear forms 



bi{ph,Vh) 
c{il^h,<Ph) 



:= -{divnh,iph)j-^ +{[fJ'h] ■n,{iph})^^, 

:= -{divfj.h,iph)j-^^ +{[tJ'h] ■n,{(ph})^i^, 

■■= {jhfi{iph),Vh).j-^^ + {rhx){^h),Lph)^^. 

It is then easy to verify that [dvfh^jbh, Vh)Th = -bi{^ih,iph) and (div/,._oA*h, 'Phjj-^ ^ 
—b{fih,Vh) for all fih G Vj^ and iph G Vj^. Furthermore, by Theorem 4.6 we have 
{Hh,'Vh'iph)j-^ = b{n h,i^h) and {fih,\'hfi'^h)j-^ = bi{fih,^Jjh). It then follows that 
we may write (6.17) in the following mixed-form: 



(6.19a) 
(6.19b) 
(6.19c) 
(6.19d) 



{vh,i'h)n +K(ih,i^h) = 
{zh,'rh)Th - K'^h,Vh) = 

-bi{zh, fh) + c{uh, tph) = 



Setting fill = Zh in (6.19a) and iph = Uh in (6.19d), we have 

{qh,Zh)rh - biizh,Uh) = 0, and - bi{zh,Uh) + c{uh,Uh) ^ 0. 

Therefore, subtracting the two equations we get {qh,Zh)Th ~ c{uh,Ufi) ~ 0. Next, 
we set Th = qh in (6.19c) and tph = Vh in (6.19b I to obtain 

{zh,qh)rh - HQh,Vh) = 0, a-nd {vh,Vh)rh + Klh^'^h) = 0. 

Adding the two equations yields ||f/i||^2(f2-) = ^{^hi Qh)Th — ~'^i'^h,Uh) < 0. There- 
fore, Vh = and c{uh,Uh) = 0. In particular, Uh vanishes on all of the boundary 
edges. Since Vh = 0, we also have 2; = by (6.19c). Setting /x/i = qh in (6.19a), 
and ij^h — Uh in (6.19b), we have 

- bi{qh,uh) 0, 
b{qh,Uh) = 0. 

Since Uh vanishes on the boundary edges, b{qh,Uh) = bj{qh,Uh). It then easily 
follows that qh = 0. Finally, we have bi{fih,Uh) = for all fih G V^!^. This in turn 
implies that 

= {nh.yuh)Th - ({M/i} • n, [m/i])^^^ = {nh,^Uh)TH V/x,i e V^. 

Therefore, Vw^jx = on all K E Th- Since [uh]\e — across all edges, we conclude 
that Uh = 0. □ 
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Remarks 6.1. 

(a) To obtain optimal order error estimates, we expect that the penalty param- 
eters must scale like r/i = 0{h~^) and 772 = 0{h^^). 

(b) The construction of DG schemes with other types of boundary conditions 
can easily be constructed by specifying the boundary data to different discrete 
differential operators. For example, if simply supported plate boundary con- 
ditions u = g and Au = q are provided, then the corresponding discrete bi- 
harmonic operator is Ah.qAh^g = divhy h,qdrvh^ h,gUh- On the other hand, 
if Cahn-Hilliard-type boundary conditions du/dn = g and dAu/dn — q are 
given, then the discrete biharmonic operator is div/i gV/idiv^ gV^u^, where 
q = qn and g = gn. 

6.3. Quasi-linear second order PDEs. 

6.3.1. The p-Laplace equation. We now extend the discrete differential frame- 
work to some non-linear elliptic problems. Although we can formulate the method 
for a very general class of quasi-linear PDEs, we shall focus our attention to a 
prototypical example, the p-Laplace equation {2 < p < 00): 

(6.20a) -div {\Vu\P^^Vu) = f in Q, 

(6.20b) u = g on dn. 



Similar to the Poisson problem, the DG method for (6.201 is obtained by simply 
replacing the grad and div operators by their respective discrete versions and adding 
a stability term to ensure that the resulting bilinear form is coercive over V/*. In 
this case, the discrete problem reads: find Uh € Vj} satisfying 

(6.21) -diVh{\W h.,gUhr^Wh,gUh) + = T^rf, 

where j)fg(-) is defined by 

(6.22) (jtjM:^'^)r. -(^i|Hr'H,b/.])£. + (^ikr'(^-5),^/.)£B 

for all iph G Vj^. Here, 771 > is a penalization parameter. Similar to the discrete 
Poisson problem, the discretization has several interpretations. By the definition 



of the discrete divergence and gradient operators, we can write (6.21 1 in the mixed 
formulation 

(6.23) V^.)r. - (Hi'^r'^'^} ' I^j), + {jhlM^V'^)n 



where q G Vj^ satisfies (6.4a I. We emphasize that the gradient appearing in the 



left-hand side of equation (6.23) is the piecewise gradient. 

Remark 6.5. In 6], Burman and Ern proposed and analyzed the following LDG 
method for the p-Laplace equation (with r ~ \ and g — Q): 

(6.24) (|V,i,oW/ir"^V/i,oW/i, V/,.,o¥'h)-^^ + (?7i|[wft]r '^['^h], VPh]) = {f,Vh)n 
for all ipfi G V^/'. Here, the authors showed the existence and uniqueness of the DG 



method (6.24 1 provided r/i = 0{h^^P). In addition, Burman and Ern showed that 
the approximate solutions converge to u strongly in LP(fl), and VhfiUh converges 
to "S/u strongly in LP{fl). Furthermore, in the one dimensional setting, numerical 
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experiments indicate a convergence rate of at least h^^^ for p e {3, 4, 5} and smooth 
exact solution. 

Clearly, method (6.24) has a similar structure to (6.23), hut they are different 
methods when p ^ 2. Indeed, since \V h.gUhY'~'^'^ h.gUh ^ , we cannot use Theo- 
rems |4.5| - |4.6| and simply write 

-{<^lYh{\V h^gUh\^^'^V h,gUh),^h).y^ = (|V/i,gW/i|^"^V/i,gUft, Vft,OV'/l)rh- 

In the following section, we show by way of numerical experiments that the DG 
method (6.23) converges with optimal order provided the exact solution is sufficiently 
smooth. 

6.3.2. Numerical experiments of the p-Laplace equation. In this subsec- 
tion we perform some numerical experiments to gauge the effectiveness of the DG 
method (6.21 ). We take the domain to be the unit square J7 = (0, 1)^ and choose 
the data / such that the exact solution is u = sin(7ra;i) sin(7ra;2) and p — 5. In all 
numerical experiments, we take the penalty parameter to be r/i = 20/hP~^ = 20//i'*. 

The resulting errors in the cases r = 1 and r = 2 are recorded in Table [T] The 
table clearly suggests that the following rates of convergence hold for smooth test 
problems: 

\\u-u^\\L2^n)=0{h^+'), 



\Vu- 



Table 1 . The errors of the computed solution and rates of conver- 



gence of the DG method (6.21 ) with solution u — sin(7ra::i) sin(7ra;2) 
on the unit square with p = 5 and t] = 20/ /i^. 





h 


\\U — UhWi.^ 


order 


|Vw - V/i«/i ^2 


order 


r = 1 


l.OOE-01 


6.48E-03 




2.33E-01 






5.00E-02 


1.Q8E-03 


2.59 


1.16E-Q1 


1.01 




2.50E-02 


1.95E-04 


2.47 


5.86E-02 


0.99 




1.25E-02 


4.99E-05 


1.97 


2.94E-02 


0.99 


r = 2 


l.OOE-01 


2.25E-03 




1.36E-02 






5.00E-02 


2.83E-04 


2.99 


3.01E-Q3 


2.18 




2.50E-02 


3.60E-05 


2.97 


7.3QE-Q4 


2.04 




1.25E-02 


4.51E-06 


3.00 


1.8QE-Q4 


2.02 



6.4. Second order linear elliptic PDEs in non-divergence form. As a fourth 
example, we consider second order elliptic PDEs written in non-divergence form. 
Namely, we consider the problem of finding a strong solution satisfying 

(6.25a) -A : D'^u ^ f in Q, 

(6.25b) u — g on dVl. 

Here, A S [C°'"(r2)]''^'^ (a € (0, 1)) is a given positive definite matrix, and A : D^u 
denotes the Frobenius inner product, i.e., A : D'^u — ^2"! j=i dx dx • 

We note that if A is sufficient smooth, e.g. if divA e L°°{^), then we may 



write the PDE (6.25a) as — div(AVu) + (divA) • Vm — f. We can then apply any 
of the standard numerical methods for convection-diffusion equations to problem 
(|6.25[). On the other hand, if A only has the regularity A G [(^o.ajdxd^ i\\en this 



argument fails and the construction of numerical methods is less obvious. As far as 
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we are aware, only two finite element methods have appeared in the literature that 
addressed the numerical approximation of problems such as (6.251. In |23j Jensen 
and Smears propose a Pi finite element method for the Hamilton-Jacobi-Bellman 
equation. To handle the lack of regularity of the coefficient matrix, the authors 
"freeze the coefficients" element-wise, and then perform the usual integration-by- 
parts technique. By modifying the framework of Barles-Sougandidis [3 , Jensen and 
Smears show that the numerical solutions converge to the exact solution strongly 



in H^. Another finite element method for problem (6.25), which is closely related 
to ours, is the one proposed by Lakkis and Pryer in [M^. Here, the authors used 
the notation of a discrete Hessian and rewrite problem (6.25) in a mixed form. 



The advantage of their approach is that the finite element spaces are simply the 
globally continuous Lagrange elements, which are simple to implement. A possible 
disadvantage of their approach is that the notion of their discrete Hessian is not 
local, and therefore writing the problem in its primal form results in a dense stiffness 
matrix. 

To formulate the DG method for the PDE using the discrete differential calculus 
framework, we again replace the continuous differential operators by the discrete 
ones. In addition, we have to project both sides of the equation onto the finite 
element space. This then leads to the following problem: find a function Uh £ 
satisfying 



(6.26) 



with ih,g{') defined by (6.3). Equivalently, the DG method is to find Uh G Vy' such 
that 



where rh G Vh satisfies (6.4a|-(6.4bl 



Remark 6.6. // the coefficient matrix A is constant, then the DG method (6.26) 
reduces to the LDG method for the PDE — div(AVu) — f with appropriate boundary 
conditions. 



Below we present some numerical test results on the DG method (6.26 1 with the 
following parameters: = (—0.5, 0.5)^, / = and 



•27) g 





if X2 


= 0, 




if xi 


= 0, 


4/3 -, 
Xi — 1 


if X2 


= 1, 


1 4/3 


if xi 


= 1 



2/3 
-1/3 1/3 



1/3 ,1/3N 
^1 ^2 
2/3 



4/3 



4/3 



It can readily be checked that the exact solution is given by u 
C^'i(ri). We note that this is a particularly challenging example since A is not 
uniformly elliptic. The resulting errors for decreasing values of h are listed in Table 
[2j and a computed solution and error is depicted in Figure [ij The table clearly 
indicates the convergence of the method, although the exact rates of convergence 
are not obvious. 



6.5. Fully nonlinear time dependent first order PDEs. As a fifth example, 
we consider Hamilton-Jacobi equations. Namely, we consider the problem of finding 
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Table 2. The errors of the computed solution and rates of con- 
vergence with ?■ = 1. 



h 


\\u - UhWl^ 


order 


|Vw - VhUhUl^ 


order 


l.OOE-01 


5.17E-03 




1.45E-01 




5.00E-02 


3.49E-03 


0.56 


9.52E-02 


0.60 


2.50E-02 


2.59E-03 


0.43 


6.50E-02 


0.55 


1.25E-02 


2.08E-03 


0.32 


4.81E-02 


0.43 



A 8.5288X10"^ 




Figure 1. Computed solution (height) and error (surface) of the 
DG method (6.26) with data (6.27) and parameters rj = 20, r — 1 
and h = 0.0125. 



the viscosity solution u € Ac C°(i}x (0,T]) for the PDE problem 

(6.28a) ut + H ( Vm) = in 17 x (0, T] , 

(6.28b) u = g on T C dQ x {0,T], 

(6.28c) u = uo on X {0}, 

where the operator _ff is a continuous and possibly nonlinear function and A is a 
function class in which the viscosity solution u is unique. We note that the following 
scheme can also be adapted for H a function of u, x, and t. 

Let B{n), BUC{n), USCifl) and LSC{n) denote, respectively, the spaces of 
bounded, bounded uniformly continuous, upper semi-continuous and lower semi- 
continuous functions on fl. We now recall the well-known existence and uniqueness 
theorem for the corresponding Cauchy problem which was first proved in 

Theorem 6.5. Let H e C(R''),, uq e BUCCR"^). Then there is exactly one 
function u G BUC{R!^ x [0,r]) for all T > such that u{x,0) — uq{x), and for 
every (fi E (R'^ x (0, oo)) and T > 0, if {xo,to) is a local maximum (resp. local 
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minimum) point of u — cj) on R'' x (0,T], then 

(f>t{xo,to)+HiV^{xo,to)) <0 
( resp. (t)t{xo,to)+H{\/<j){xo,to)) > 0). 

Definition 6.1. The function u whose existence and uniqueness is guaranteed by 



Theorem 6.5 is called the viscosity solution of the Cauchy version of (6.281 



Recently, a nonstandard LDG method was proposed by Yan and Osher in 



for approximating the viscosity solution of the Hamilton- Jacobi problem (6.28). 
The main idea of [34] is to approximate the "left" and "right" side derivatives of 
the viscosity solution and to judiciously combine them through a monotone and 
consistent numerical Hamiltonian (cf. [33]) such as the Lax-Friedrichs numerical 
Hamiltonian 

(6.29) Hiq-,q+) := H (^^^) - • {q+ q-) , 

where f3 E is an undetermined nonnegative vector chosen to enforce the mono- 
tonicity property of H, or the Godunov numerical Hamiltonian 

(6.30) H(q^ ,q^) -.^ ext - +\---ext - +\H(q), 
where 

_ \mma<v<b, if a < &, 
ext„g/(Q M :— < -r 1 

|^maxb<^,<a, it a > &, 

for I{a,b) := [min(a, 6), max(a, 6)] . 

Remark 6.7. For r — and (3 — \ on a uniform rectangular grid, the second 



term in (6.29) is equivalent to a second order finite difference approximation for 



the negative Laplacian operator scaled by h. Thus, the second term in (6.29) is 
called a numerical viscosity, and the method is a direct realization of the vanishing 
viscosity method from PDE theory. However, for high order elements and variable 
coefficient vector f3, while we do not exactly recover a scaled Laplacian operator, 
we do recover some of the stabilizing properties from adding a second-order-like 
perturbation. 

With the correct choice of discrete derivatives, we can rewrite the nonstandard 
LDG method of Yan and Osher [34] as follows: For each time step n = 1,2,3,... 
with u^J := V^uq, find uJJ G using the recursive relation 



ul-' + ^tV^H{V-ul-\^lul 



In addition to the above Euler time-stepping method with H given by the Lax- 



Friedrichs numerical Hamiltonian defined in (6.29), in [34] Yan and Osher also im- 
plemented the explicit third-order TVD Runge-Kutta time-stepping method given 
in [3D]. Tests include one and two dimensional problems for r > 0. 



Remark 6.8. For r — on a Cartesian grid, the method reduces to the convergent 
FD method of Crandall and Lions proposed and analyzed in |12j . 
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6.6. Fully nonlinear second order PDEs. As the last example, we consider 
fully nonlinear second order elliptic PDEs. Namely, we consider the problem of 
finding the viscosity solution u E A C C'^(51) for the PDE problem 



(6.31a) 
(6.31b) 



F[u] -.^ F (D^ 



in fl, 
on dil, 



where the operator F can be nonlinear in all arguments and A is a, function class in 
which the viscosity solution is unique. Throughout this subsection, we assume that 



F[v] is elliptic for all v G A, F satisfies a comparison principle and problem (6.31) 
has a unique viscosity solution u € A. The definitions for the above terms are given 



below. We again use the same function and space notations from Section [63] We 
also note that the following scheme can also be adapted for F a function of u and 
Vu. 



For ease of presentation, we write (6.31 ) as 
.32) F{D^u,x) 







in il, 



where we have used the convention of writing the boundary condition as a discon- 
tinuity of the PDE (cf. [31 p.274]). Also for ease of presentation, we assume that 
F is a continuous function and refer the reader to O [T71 [TH] for the case when F 
is only a bounded function. 

The following three definitions are standard and can be found in [HI [7l [3] . 



Definition 6.2. Equation (6.32) is said to be elliptic if for all x £ there holds 

(6.33) F{A, x) < F{B, x) VA, B e S"^""^, A>B, 

where A> B means that A — B is a nonnegative definite matrix, andS'^^'^ denotes 
the set of real symmetric d x d matrices. 

We note that when F is differentiable, ellipticity can also be defined by requiring 



that the matrix 



OF 
dD^u 



is negative semi-definite (cf. |21j p. 441]). 



Definition 6.3. A function u G C^iVl) is called a viscosity subsolution (resp. 
supersolution) of (6.32) if, for all if £ C^(il), if u — ip (resp. u — ip) has a local 
maximum (resp. minimum) at xq G fl, then we have 

F{D^^{xo),xo) < 

(resp. F{D'^(p{xo),Xo) > 0). The function u is said to be a viscosity solution of 



(6.32) if it is simultaneously a viscosity subsolution and a viscosity supersolution of 



(6.31) 



Definition 6.4. Problem (6.32) is said to satisfy a comparison principle if the 

following statement holds. For any upper semi- continuous function u and lower 
semi-continuous function v on fl, if u is a viscosity subsolution and v is a viscosity 
.supersolution of (6.32), then u <v onfl. 



Inspired by the work of Yan and Osher [51], the first and second authors of 
this paper recently proposed in |17| a nonstandard LDG method for approximat- 
ing the viscosity solution of the fully nonlinear second order problem (6.31) in 
one-dimension. The main idea of [17] is to use all four of the various "sided" ap- 



proximations for the second order derivative, (3.13a), of the viscosity solution and 



to judiciously combine them through a g-monotone (generalized monotone) and 
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consistent numerical operator such as the following Lax-Friedrichs-likc numerical 
operator that has been adopted from [TTI for an arbitrary dimensional problem: 

(6.34) F(P— ,p-+,P+-,P++,0 ^^^^ ,e 

+ A : {P— - P + - P+ + P++), 

where A £ Jid'^d jg ^ nonnegative constant matrix that is chosen to enforce the 
g-monotonicity property of F. The consistency of F is defined by fulfilling the 
following property: 

F{P,p,p,p,x) = F{P,x) yp e'R'^'"^,x en; 

and the g-monotonicity requires that F is monotone increasing in its first and 
fourth arguments (i.e., P , P^"*") and monotone decreasing in its second and 
third arguments (i.e., P ^, P^ ). 

Below we give a reformulation of the nonstandard LDG method of |17j using our 
DG finite element differential calculus machinery. To this end, we simply replace 
the continuous differential operators by multiple copies of discrete ones. Due to 
the lack of integration by parts caused by the nonlinearity, we have to project the 
equation onto the DG finite element space. This then leads to the following scheme 
of finding u/j G such that 



(6.35) P/? F [D--gU,,, D-+u„, D+'^u^, D++u„, xj - 0, 

where , ^ht> -^/lo' ^'-'^^ sided numerical Hessians (with the 



h,g ' h,g ' h,g ' h,g 

prescribed boundary data g) defined in Sections] It can be shown that (6.35) is 
indeed equivalent to the LDG scheme of [17j in one-dimension. 

Remarks 6.2. 



(a) When r — and d — 1, scheme (6.35) reduces to the FD method given in 



|18j . which was proved to be convergent. 
(b) For r = and A — Idxd on a rectangular grid, the second term on the 



right hand side in (6.34) is equivalent to a second order finite difference 



approximation for the biharmonic operator scaled by h^ . Thus, the second 



term in (6.34) is called a numerical moment, and the method is a direct 
realization of the vanishing moment method proposed in [191120) . However, 
for high order elements and variable coefficient matrices A, while we do 
not exactly recover a scaled biharmonic operator, we do recover some of the 
stabilizing properties from adding a fourth- order-like perturbation. 

(c) Numerical tests of |17j show the above discretization can eliminate many, 
and in some cases all, of the numerical artifacts that plague standard dis- 
cretizations for fully nonlinear second order PDEs ( cf. | 16| and the refer- 
ences therein). 

(d) Fully discrete schemes of [T7| for parabolic fully nonlinear second order 
PDEs can also be recast using the DG finite element differential calculus 
machinery. 



To solve the algebraic system (6.35 1, a nonlinear solver must be used. Numerical 
tests of [17] show that when the initial guess for Uh is not too close to a non- 
viscosity solution of the PDE problem, a Newton-based solver performs well as 
long as the solution is not on the boundary of the admissible set A. However, in 
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the degenerate case, a split solver based on the DWDG discretization for the Poisson 



equation from Section 6.1.2 and the Lax-Friedrichs-like operators in (6.34) appear 



to be better suited. This new solver for (6.35) is given below in Algorithm 6.1 



We note that this solver appears to work well even for some cases when the initial 
guess for is not in A. Thus, the solver uses key tools from the discretization to 
address the issue of conditional uniqueness of viscosity solutions. 



Let A^gt; and A^^w denote the diagonal ma- 



Algorithm 6.1. Pick u^^^ e Vj^ 
trices formed by the diagonals of D^^v and D^^v, respectively, and A^^w and 
A^"^?; denote the diagonal matrices formed by the diagonals of DJ^~v and D'^'^v, 

respectively, for all v e 1//' . Let XI denote the diagonal matrix formed by the vector 
A G R''. 

For n = 1,2,3,..., 

Step 1: For i = 1,2, . . . ,d, set 



\GP := f[{{D-+ A-+)/2 + {D+-^ A+; 
+ A 



)/2)^ 



A(«)I,: 



(n-l) 

h 



A, u 



A++4""'^ -2A(")I 



for a fixed constant A > Q such that G'-"-' is monotone decreasing with 
respect to A*-"^ . 
Step 2: Solve for A^") G such that 

([Gt\ 



Th 



for all 4>i G Vj}, i = \,2,...,d. 
Step 3: Solve for u-^^ G such that 



2 (^/i.s"/! 



(") V7 + 



h,0 



(n) 



d 

E 



(A 



Th 



for all Vh G T^/* . 

Note that we are solving the discretization that results from the choice A — Aldxd 



in (6.34) 



The above solver is a fixed point method for the diagonal of the Hessian approx- 
imation formed by {D^~^ + D^~)/2. We can see that the nonlinear equation in 
Step 2 is entirely monotone and local in its nonlinear components when A is suffi- 
ciently large. Step 3 is well-defined due to the well-posedness of the DWDG method 
for Poisson's equation. Lastly, as written, we can see that the numerical moment 
inspired by the vanishing moment methodology serves as the motivation for using 
the diagonal of the Hessian approximation as the fixed-point parameter. Numerical 
tests indicate that the solver destabilizes numerical artifacts even when they exist 
(cf. [TT]). Thus, the solver directly addresses the issue of conditional uniqueness by 
enforcing the preservation of monotonicity in the Hessian approximation at each 
iteration. 
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